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Abstract 

Analyzing diverse seismic catalogs, we have determined that the probability 
densities of the earthquake recurrence times for different spatial areas and 
magnitude ranges can be described by a unique universal distribution if the 
time is rescaled with the mean rate of occurrence. The shape of the distri- 
bution shows the existence of clustering beyond the short-term regime, and 
scaling reveals the self-similarity of the clustering structure. This holds from 
worldwide to local scales and for quite different tectonic environments. After- 
shock sequences also follow this universal recurrence-time distribution if the 
rescaling is undertaken with the time-varying rate. 
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Although earthquakes are a phenomenon of great complexity, certain simple general laws 
govern the statistics of their occurrence (1-5); however, no such a unified description has 
yet been established for the time interval between successive events (6-9). On the contrary, 
the rich variability intrinsic to earthquakes has promoted that all possibilities have been 
proposed, from totally random occurrence to the periodic ticking of great quakes. The most 
extended view is that of two separated processes, one for mainshocks, which ought to follow a 
Poisson distribution (10) (or not (6-9)), and an independent process to generate aftershocks. 
Consequently, the "standard practice" for this approach consists first of delimiting the spatial 
area to be studied, on the basis of its tectonic characteristics, and then of a careful (and not 
so standard) identification of aftershocks, in order to separate them from the main sequence. 

Here we take an alternative perspective, complementary to the previous reductionist 
view. We try to look at the system as a whole, irrespective of tectonic features and placing 
all the events on the same footing, whether these would be classified as mainshocks or 
aftershocks (11,12). This follows one of the key guidelines of complexity philosophy, which 
is to find descriptions on a general level; the existence of general laws fulfilled by all the 
earthquakes will unveil a degree of unity in an extremely complex phenomenon (13). 

We analyze a global catalog, the PDE from the NEIC ( 14 ), as well as several local cata- 
logs: that of the SCSN (Southern California) (15), the JUNEC (Japan) (16), the Bulletins 
of the IGN (the Iberian Peninsula and the North of Africa) (17), and the BGS catalog (the 
British Islands and the North Sea) (18). Catalogs generally characterize each earthquake by 
three main quantities: time of occurrence, magnitude, and a vector of spatial coordinates 
for the hypocenter. These are also the variables that we focus on. 

Without concerning ourselves with the tectonic properties, following Bak et al. (11,12), 
we consider spatial areas delimited by a window of L degrees in longitude and L degrees in 
latitude (this corresponds to a square region if these angles are translated onto a rectangular 
coordinate system (19)). For each one of these regions, only events with magnitude M above 
a certain threshold M c are considered (the threshold should be larger than the minimum 
magnitude for which the catalog is considered complete). In this way, we transform a 
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time process in four dimensions (spatial coordinates and magnitude) into a simple process 
on a line for which events occur at times t iy with % = 1,2... TV, and therefore, the time 
between successive events can be obtained as Tj = t i+ i — U,i = 1,2 ... N — 1. These are 
the recurrence times in a given L 2 — region for events above M c , which can also be referred 
to as inter-occurrence or waiting times. Note that, with this transformation, we have lost 
the structure in space and in the magnitude scale; nonetheless, the change in the process 
properties with the variation of L and M c will allow us to recover some of this information. 

Due to the multiple time scales involved (from seconds or minutes to many years), the 
probability density of the recurrence time must be calculated with care. We could work 
with the logarithm of r, but an equivalent and more direct possibility is to define the bins 
over which the probability density is calculated exponentially growing as c n , with c > 1 
and n labeling consecutive bins (we usually take c = 2.5). This ensures that we have the 
appropriate bin size for each time scale. We then count the number of pairs of consecutive 
events separated by a time whose value lies into a given bin, and divide by the total number 
of pair of events (number of events minus one) and by the size of the bin, to attain the 
estimation of the probability density D xy {r) over that bin, where xy denotes the spatial 
coordinates of the region. (D xy also depends on L and M c , but for the sake of simplicity 
in the notation, we obviate this dependence.) Moreover, due to the incompleteness of the 
catalogs in the short-time scale, we will not display in the plots time intervals that are 
smaller than 2 minutes. 

The entire Earth has been analyzed by this method. Figure 1A shows the results for 
D xy (r) for worldwide earthquakes in the NEIC-PDE catalog for the 1973-2002 period, using 
L from 180° to 2.8° (about 300 km) and M c from 5 to 6.5. Note the variation of the 
recurrence time across several orders of magnitude. Figure IB shows the rescaling of all 
the distributions with the mean rate R in the region, defined as the total number of events 
divided by the total time interval over which these events span. The perfect data collapse 
implies that we can write 
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where R xy stresses that the rate refers to the (x, y)-region (of size L 2 and with M > M c ). 
The scaling function / can be well fit by a generalized gamma distribution, 

f(9) = C^eM-d 5 /B), 

with parameters 7 = 0.67±0.05, 5 = 0.98±0.05, B = 1.58±0.15, and C = 0.50±0.10, which 
has a coefficient of variation CV ~ 1.2. In fact, the value of delta can be approximated to 
one, which corresponds to the standard gamma distribution; we therefore essentially have 
a power law with exponent about —0.33, up to the largest values, where the exponential 
factor comes into play. 

The scaling function fits the rescaled distributions surprisingly well for intermediate and 
large values of the recurrence time, about r > 0.01/it^, (this usually contains from 90% 
to 95% of probability). The deviations are considerable for small values of r; although 
the statistic is low in this case (few events in the small bins being considered), for certain 
regions there is a clear tendency for the distribution to exceed the value given by the scaling 
function, that is, there is an excess of very short recurrence times, in the form of another 
power law but with the exponent close to (minus) one. This occurs when the rate in the 
region is not stationary, due to the sudden increase and slow decay of the activity provoked 
by aftershock sequences. Furthermore, these increments become more apparent when the 
size of the region decreases, in such a way that, for L < 11°, not all the regions in the world 
verify the scaling law; this can be solved in some cases by rescaling with the mean rate in 
the region calculated, not over the whole time-span of the catalog, but only over the period 
for which the rate is stationary (no activity peaks). Nevertheless, the aftershocks can be so 
important for certain particular regions that a stationary period may not exist. 

The same analysis is performed on local catalogs and identical results hold, as Fig. 1C 
illustrates. For Southern California for the 1984 — 2001 period, a number of small regions 
with stationary activity are shown; for larger regions, the time window must be reduced to 
1988 — 1991 or to 1995 — 1998, for instance, in order to find stationariness. For Japan, we 
also analyze large regions for the 1995 — 1998 period, for the Iberian Peninsula the period 
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is 1993 — 1997, and for the British Islands, 1991 — 2001. The magnitude thresholds range 
from 2 to 4, and L from 30° to 0.16° (approximately 3300 to 17 km). 

The shape of the distribution D xy (r) indicates that the memory of the last earthquake 
is conserved up to the largest times, with the probability of a subsequent event being maxi- 
mum immediately after the last shock, and slowly decreasing with time. This constitutes a 
clustering effect (20), in which earthquakes attract each other, and has as a counter-intuitive 
consequence the fact that the longer it has been since the last earthquake, the longer the 
expected time will be till the next (21). On the other hand, the scaling of D xy (r) under 
changes in M c , L, and the region coordinates implies that the clustering structure is self- 
similar over different regions and magnitude ranges. The robustness of the distribution under 
such changes is therefore noteworthy. It is also remarkable that if the region is kept fixed 
and only M c varies, the scaling with the rate R xy can be substituted by the factor 10~ bMc , 
where b refers to the b— value of the frequency-magnitude relation (1,2) in that particular 
region (11,12); despite the regional variability of b, the universality of the scaling function 
/ remains valid. 

At the outset of this exposition, we suggested that our results were valid for all events, 
including aftershocks; however, the aftershocks can break the scaling of the distribution up 
to very large time values. We can then say that the universal distribution does not describe 
the short-time intervals over which aftershocks replicate. This may appear to be correct, 
but can be turn out to be mistaken in the following way: for aftershock sequences, the mean 
rate is not stationary, but changes with time; therefore, we should rescale the recurrence 
times using the "instantaneous" rate r xy (t). For many sequences, r xy (t) is found to decay 
following the Omori law (3), 

r xy{t) — ^ p 

where t is the time elapsed since the mainshock (and the parameters depend on x, y, L, 
and M c ). Figure 2 A shows precisely this for several important earthquakes in Southern 
California (15). The corresponding recurrence-time distributions, calculated over the period 
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for which the Omori law is fulfilled, are displayed in Fig. 2B; the results after rescaling with 
r xy (t) appear in Fig. 2C, again in agreement with the universal distribution. The scaling is 
outstanding, taking into account that the p— values spread from 0.9 to 1.35. Note also that 
this scaling implies the existence of a secondary clustering inside the primary clustering of 
the aftershock sequence, and therefore the process is not a nonhomogeneous Poisson process. 

The present characterization of the stochastic spatio-temporal occurrence of earthquakes 
by means of a unique law would indicate the existence of universal mechanisms in the 
earthquake-generation process (22), the understanding of which, however, is still far be- 
yond us. Nevertheless, the context of self-organized critical phenomena (11-13) provides a 
coherent framework at this stage. These findings can also be relevant to continuous-time 
random- walk models of seismicity (23,24), time-dependent hazard, and forecasting in general 
(25-28). 
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Figure Captions 

Fig. 1. Recurrence-time distributions without and with rescaling. (A) Probability 
densities from the NEIC-PDE worldwide catalog for several regions, L, and M c . For L > 
45°, all the regions with more than 500 events are shown, whereas for L < 22.5°, only 
representative regions with moderate aftershock activity are displayed. The vector (k x ,k y ) 
labels the different regions, for which the coordinates of the center can be obtained as 
x = x min + (k x + 0.5)L, y = y min + (k y + 0.5)L, with x min = -180°, y min = -90°. (The 
360° x 180° region, which covers the whole planet, has been included for completeness.) (B) 
Previous data, after rescaling, with a fit of the scaling function /. (C) Rescaled distributions 
from local catalogs. SC88, SC95, SC84 refer to Southern California for the years 1988-1991, 
1995 — 1998, and 1984 — 2001. To obtain region coordinates, use the previous formula with 
(x m in,ymin) = ("124°, 29°), (-123°, 30°), (125°, 25°), (-20°, 30°), and (-10°, 45°) for SC88- 
95, SC84, Japan, the Iberian Peninsula, and the British Islands respectively. The function 
displayed is the fit obtained from the NEIC-PDE catalog. 

Fig. 2. Analysis of aftershock sequences. (A) Decay of the rate after a mainshock and 
illustration of the Omori law for the following earthquakes in Southern California: Chalfant 
Valley (July 21, 1986, M = 5.9), Landers (June 28, 1992, M = 7.3), Northridge (Jan. 17, 
1994, M = 6.7), and Hector Mine (Oct. 16, 1999, M = 7.1). Regions of diverse size L are 
considered, all of these including the mainshock. Some curves are shifted for the sake of 
clarity. (B) Distributions of recurrence times for the previous sequences. (C) Distributions 
ip of the dimensionless time r xy (t)r, in total agreement with the universal scaling function 
/■ 
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